/*******************************************************************
 *   Mandelbrot set computation
 *
 *   CSIS 248 
 *   A Varghese
 *   
 *   For more information on the Mandelbrot set,
 *    see: http://en.wikipedia.org/wiki/Mandelbrot_set
 ******************************************************************/

#include <stdio.h>
#include <math.h>
#include <sys/stat.h>
#include <fcntl.h>

typedef
struct complex{
   double re;
   double im;
} Complex;

int main(){
   /*
    * Lower left and upper right corners of complex plane
    */
   const Complex ll = {-2.5, -1.0};
   const Complex ur = { 1.0,  1.0};

   /*
    * Resolution: Nx by Ny pixels in complex plane
    */
   const int Nx = 40;
   const int Ny = 40;

   const double dx = (ur.re - ll.re)/(double)(Nx-1);
   const double dy = (ur.im - ll.im)/(double)(Ny-1);


   Complex z[Nx][Ny];
   Complex c[Nx][Ny];
   unsigned char pixel[Nx][Ny];

   /*
    * Number of iterations
    */
   const int number_of_iterations = 100000;
   /*
    * Once the real or imaginary parts go beyond the max
    *  values, they have "escaped"
    */
   const double Max = 4.0;

   /*
    * Initialize z and c
    */
   printf("Initializing ... \n");
   int i=0;
   for (i=0; i<Nx; i++) {
      int j=0;
      for (j=0; j<Ny; j++) {
        z[i][j].re = z[i][j].im = 0.;

        c[i][j].re = ll.re + i * dx;
        c[i][j].im = ll.im + j * dy;

	pixel[i][j] = 1; // close to black
    }
  }

   /*
    * Iterate for number of iterations
    */
   printf("Starting iterations ... \n");
   int k=0;
   for(k=0; k<number_of_iterations; k++){
      /*
       * Figure out the color for this iteration
       */
      const unsigned char color = (unsigned char) ((k * 255) / number_of_iterations);

      for (i=0; i<Nx; i++) {
         int j=0;
         for (j=0; j<Ny; j++) {
	    /*
             * if pixel value is 0, no need to compute anything.
             */
            Complex t = z[i][j];
            /*
             * check magnitude
             */
            if (t.re * t.re + t.im*t.im < Max) {
		/*
		 * Compute z_n+1 = z_n * z_n + c
		 */
		Complex znp1;
		znp1.re = t.re * t.re - t.im * t.im + c[i][j].re;
		znp1.im = 2.0 * t.re * t.im         + c[i][j].im;

		z[i][j] = znp1;

		pixel[i][j] = color;
      
	    } else {
		pixel[i][j] = 0; // black
	    }
	 }
      }
   }
   printf(" ... done!\n");

   /*
    * Dump output to screen
    * - Note that we write it out "rows first"
    *   starting with the largest y value (imaginary component)
    */
   int j=Ny-1;
   for (  ; j>=0; j--){
      for (i=0; i<Nx; i++)
         printf("%03u ", pixel[i][j]);
      printf("\n");
   }

   /*
    * Dump binary output to file
    */
   const char filename[] = "mandelbrot.pix";
   int out = open( filename, O_RDWR | O_CREAT, S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH );
   if (out < 0){
      printf("Not able to open file %s for writing\n", filename );
      return -1;
   }
   for (j=Ny-1  ; j>=0; j--){
      for (i=0; i<Nx; i++)
	write( out, &(pixel[i][j]), 1 );
   }
   close( out );
}
